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Abstract 

Based on a previously developed recursive approach for calculating the short-time expansion of 
the propagator for systems with time-independent potentials and its time-dependent generaliza- 
tion for simple single-particle systems, in this paper we present a full extension of this formalism 
to a general quantum system with many degrees of freedom in a time-dependent potential. Fur- 
thermore, we also present a recursive approach for the velocity-independent part of the effective 
potential, which is necessary for calculating diagonal amplitudes and partition functions, as well 
as an extension from the imaginary-time formalism to the real-time one, which enables to study 
the dynamical properties of quantum systems. The recursive approach developed here allows an 
analytic derivation of the short-time expansion to orders that have not been accessible before, 
using the implemented SPEEDUP symbolic calculation code. The analytically derived results are 
extensively numerically verified by treating several models in both imaginary and real time. 
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I. INTRODUCTION 



Although a large number of physical systems admit studies of their basic properties using 
different types of time- independent formalisms, in many important cases one has to explicitly 
take into account the time dependence and to apply the appropriate approach, i.e. one of 
available time-dependent analytical or numerical methods. This has to be done even for 
systems with time-independent potentials if their dynamical properties and time evolution 
is studied. However, for describing systems in genuinely time-dependent external potentials 
or for rotating systems, using of such approaches is always necessary. This applies equally 
to classical and quantum systems, and various methods have been developed to address 
relevant physical problems. For quantum systems a number of general methods is available, 
ranging from time-dependent perturbation theory and variational perturbation theory, to 
specialized approaches such as the Density Matrix Renormalization Group (DMRG) Q, Q|, 
the Density Functional Theory (DFT) [3|, |4j, and the Density Matrix Functional Theory 
(DMFT) [5|. 

In addition to these generic schemes, several specific numerical methods have been de- 
veloped for enabling a quantitative description of quantum systems that have attracted a 
significant research mterest. This includes studies of ultra-cold quantum gases and 



optical lattices 



16 



whose comprehensive and hig. 



lly tunable features make them an 



important example of Feynman's quantum simulator 



21 



201 ] . In such numerical approaches 



-|24j usually a second-order algorithm in the propagation time is used, which is basically 



the same as in the time-independent case. However, also higher order schemes have been 



derived, including fourth 



25 



- |28| and higher-order expansions 29|-|32|. 



The main goal of this paper is to develop a formalism which enables a systematic improve- 
ment in the convergence order of numerical algorithms for general time-dependent quantum 
systems. In our previous paper 33], we have extended the earlier established approach 



34J-|371] for obtaining a high-order short-time expansion of transition amplitudes for time- 



independent potentials to the important time-dependent case. In that paper, we have first 
calculated a short-time expansion of a generic transition amplitude using the path integral 
formalism to fourth order in the propagation time, which has then served the important 
purpose: to introduce the proper functional form of the ansatz for the ideal discretized ef- 
fective potential, as well as its short-time double series expansion in both the propagation 



time and the discretized velocity. Using both the forward and backward Schrodinger equa- 
tion for transition amplitudes, we have then derived the appropriate equation for the ideal 
effective potential, which represents an important generalization of the equation derived in 
Ref. [37J for the time-independent case. Using the double series expansion ansatz for the 
effective potential, we have been able to set up and solve recursive relations for the effective 
potential for one-dimensional systems, and to numerically verify that higher-order analytic 
approximative results for transition amplitudes obtained in this way, indeed, give the correct 
convergence order for a number of models. 

In thispaper we further develop and generalize the time-dependent approach introduced 



in Ref. 



33]. First, in Sec. HI] we briefly review the time- dependent effective action approach 



in the path- integral formalism of quantum mechanics |38H4l| . as well as the main results of 
our previous paper [33]. Then we extend the recursive approach for calculating the short- 
time expansion of the effective potential to quantum systems with many degrees of freedom 
in Sec. IIHI In this section, the generalized approach is also numerically verified for the 
case of a simple time-dependent multi-component harmonic oscillator system, and possible 
relevant physical applications in the realm of ultracold quantum gases are briefly indicated. 
In Sec. IIVI we present another important extension of the time-dependent formalism, and 
set up a specific recursive relation for calculating diagonal transition amplitudes, which 
are necessary for a numeric high-precision calculation of partition functions. In this case, 
a simplified set of recursive relations is obtained, which is numerically verified for several 
model potentials. Finally, Sec. |V] illustrates how the developed imaginary-time formalism 
can be transformed into a real-time one, and its applicability is numerically demonstrated 
by treating several models. We also show how the real-time evolution can be described using 
the time-dependent effective action approach, and how the associated numerical errors can 
be assessed and controlled in typical applications. 



II. EFFECTIVE ACTION APPROACH FOR SYSTEMS WITH TIME-DEPEN- 
DENT POTENTIALS 



In this section we give a brief overview of the time-dependent effective action approach es- 
tablished in Ref. 33] . This formalism considers a non-relativistic quantum multi-component 

3 



system in d spatial dimensions with a Hamiltonian of the form 



P 



£(p,q,0 = E^+ y (M> (2.1) 



i=l 



2M 



(0 



where P denotes the number of particles in the system and the P x d dimensional vectors 
q and p contain positions and momenta of all particles, while the parenthetic subscript (i) 
denotes the particle number. For such a system we consider the calculation of the transition 
amplitudes 

A(a,t a ;b,t b ) = (b,t b \U(t a ^t b )\3L,t a ), (2.2) 

where the vectors a and b describe the positions of all particles at the initial and final time 
t a and tb, respectively. The above transition amplitude is a coordinate-space matrix element 
of the evolution operator, which describes the propagation of the system (12.11) from t a to t b , 
and is defined by the time-ordered exponential 

h 



U{t a ^t b )=Texp\-l I dtH{p,q,t)}. (2.3) 



t a 

In the path-integral formalism, the transition amplitude can be expressed by a coordinate- 
space path integral 

r q{t b )=b 
«/q(ta)=a 

where the integration is defined over all possible trajectories q(t). This usually involves the 
discretization of the trajectories, which is usually performed by dividing the time evolution 
from t a to t b into iV equal time steps. In the above equation, 5* denotes the action for a 
given trajectory q(t): 



.,4(a./-„:b. /,,) = / Pq(/) <>xp <{ ^S[q\ } . (2.1) 



1 



S[q] = jf j-q 2 (t)-y(q(t),t)j. (2.5) 

The common step at this point is to switch to the imaginary-time formalism, which is 
usually applied in numerical simulations [42], due to problems which may arise from the 
oscillatory nature of the integrand in the real-time approach. We will do so as well in the 
next two sections, but in Sec. |V] we will switch the developed formalism back to the real 
time and demonstrate how it can be used for studying the dynamics and real-time evolution 



of quantum systems. After Wick rotation to the imaginary time, the transition amplitude 
in the path-integral formalism is expressed as 

"i(t b )=b 



A(a, t a ;b,t b 



Z>q(t) e~^M 



(2.6) 



'q(*a)= a 

where the Minkowski action is now replaced by its imaginary-time counterpart, the Euclidean 
action 



5 E [q] 



to 



dt^-ci 2 (t) + V(q(t),t) 



(2.7) 



which is actually the energy functional for the system. 



As was shown previously for time- independent potentials 371] . as well as for the time- 
dependent ones in Ref. 33], the exact imaginary-time transition amplitudes can be expressed 
in the form 



A(a, t a ;b,t b ) 



1 



-S*(x,x;e,r) 



(2.8) 



(2tt £ )™/ 2 

where S* stands for the ideal discretized action and depends on the coordinate mid-point 
x = (a + b)/2, the discretized velocity x = (b — a)/2, the time interval e = t b — t a , and the 
time mid-point r = (t a + t b )/2. Note that we have used the convention h — 1, and we have 
restricted ourselves to particles with unitary masses mu) = 1. The ideal discretized action 



further reads 



33 
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S^x, x; e, r) = -x 2 + eWhc, x; e, r) 

e 



(2.9) 



where is the ideal discretized effective potential, which also depends on the time mid- 
point t due to the explicit time dependence of the potential V. This represents the major 
difference in the formalism compared to the previously developed one in Ref. [37j for the 
time-independent case. 

In order to determine a partial differential equation for the ideal effective potential W, we 



have derived in Sec. 4 of our previous paper 



33| the forward and the backward Schrodinger 



equation for time derivatives of the transition amplitude with respect to the initial and final 
time, and have found that they obey 



d tb A(&, t a ; b, t b ) = -H b A(a, t a ; b, t b ) 
d ta A(a, t a ; b, t b ) = H a A(a, t a ; b, t b ) , 



(2.10) 
(2.11) 



where Hj, stands for the coordinate-space Hamilton operator Hj, = H{— id^,, b, t b ), in which 
momentum and position operators are replaced by their coordinate-space representations 
at b (and similarly for H a ). When we insert the ansatz (12. 8p together with (12. 9p into the 
equations (12.101) and (12. lip , we obtain the corresponding partial differential equation for 
the effective potential in the form 

W+x-dW+ed e W--ed 2 W--ed 2 W+-e 2 {dW) 2 +-e 2 (dW) 2 = -(V++V-), (2.12) 
8 8 8 8 2 

where V± = V(x ± x, r ± |), i.e. V_ = V(a, t a ) and V + = V(b, U). As expected, the form of 



this equation is the same as Eq. (29) from Ref . [37| , and if the potential V does not depend 
on time, we immediately obtain the previously derived result. 

The above equation is the starting point for calculating the effective potential. Naturally, 
it can be solved analytically only for exactly solvable models. However, if the propagation 
time is short, we can perform a short-time expansion of the effective potential and set up 
appropriate equations for the coefficients in this expansion. In our previous paper {33} this 
was done for one-dimensional systems, and it was shown that we have to use a double 
expansion of the effective potential in both e and x. The reason for this is the fact that the 
propagation in the imaginary time is equivalent to diffusion, and therefore, on the average, 
we expect the diffusion relation x 2 ~ e to hold. This allows us to effectively couple the two 
expansion parameters and to establish a unique counting of powers in e for all terms in the 
expansion for W . However, if the diffusion relation is not applicable, the two expansions 
can be considered as independent, and the whole approach can still be applied, with an 
independent counting of powers in e and in x. 

In the next section we will analytically derive a systematic short-time expansion of the ef- 
fective potential W for quantum many-body systems, which yields a significant improvement 
in the convergence of numerically calculated transition amplitudes and partition functions 
for systems in time-dependent potentials. As we see from Eq. (12.91) . if the effective potential 
W is calculated to order we get the effective action correct to order e p . If we take 

into account the normalization factor in the expression ( 12.81) . the corresponding error in the 
calculation of the short-time amplitude is given by 

A p {a,t a ;b,t b ) = A^t a -h,t b ) + 0(e p+1 ~ pd l 2 ), (2.13) 

where subscript p denotes that we use the effective action of that order in e. Therefore, 
we see that the analytical calculation of higher-order effective actions is beneficial, since 



it provides an analytic approximation for transition amplitudes which yield high-precision 
results of the desired order in numerical calculations. 



III. MULTI-COMPONENT SYSTEMS 



Now we focus on the development of the recursive formalism for calculating the effective 
potential for the case of a general mult i- component non-relativistic quantum system of P 
particles in d dimensions by extending the one-dimensional calculations of Sec. 5 in our 
previous paper 33]. Note that such a formalism is needed even for studies of single-particle 



systems in two or three spatial dimensions, which have more than one degree of freedom. For 
example, such time-independent many-body effective actions of level p = 21 have already 
been used for a numerical study of fast-rotating Bose-Einstein condensates 45], as well 



as a high-precision calculation of the energy spectra and eigenfunctions of several two- 
dimensional models [46,0]. The presented extension of the many-body formalism will allow 
studies of such systems in external time- dependent potentials, as well as the investigation of 



the formation and evolution of vortices 



14 



151 ] and other dynamical phenomena. We also 



plan to study collective oscillation modes of Bose-Einstein condensates with a parametrically 



modulated interaction 
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To develop the time-dependent many-body formalism, we solve the partial differential 
equation ( 12.121) for the effective potential W by using a multi-dimensional many-particle 
generalization of the double power expansion used in Eq. (33) of our previous paper |33|] for 
one-component systems, which has the form 



oo in 

W(x, x;e,T) = £E { W ™A^ *; r) e m ~ k + W m+1/2jfc (x, x; r) s m - k } . 



(3.1) 



m=0 k=0 

Here we have introduced the contractions 



W m , fc (x,x; r) 



— li,...t'2k ( . \ 



(3.2) 



WWi/2,fc(x,x; t) 



f. ...f. *i,-»a*+i 
x *l x *2fc+i SrH-l/2,fc 



x; r) 



(3.3) 



in such a way that they correspond to the case of the time-independent potential 37J • In the 
above relations we assume the Einstein convention that summation over repeated indices is 
performed. Introducing such contractions of tensorial coefficients c significantly simplifies 



in this case the analytic derivation and also provides a key ingredient for implementing the 
many-body recursion relations in symbolic calculations using e.g. Mathematica software 



package 50]. Otherwise, the task to explicitly symmetrize the coefficients would amount to 
a complexity of the algorithm which scales with the number of possible permutations (Pxd)\ 
and which would not be feasible even for a very moderate number of particles. Using scalar 
quantities, which are obtained by contracting the coefficients with the discretized velocity 
x, efficiently solves this problem. 

The multiple-component version of the expression from the right-hand side of Eq. (I2.12p 
has the form 

YYI n(m,AQ e "-* 2k (1-U(m,k))e^ , ± . a) »*i\ % k) (3 4) 

^ ^l(2fc)!(m-fc)!2— k[ ' + (2k + l)\(m-k)\2^ [ °> J ' 

(m—k) 

where V represents (m — fc)-th partial derivative of the potential with respect to the time 
t. After inserting the above expression into the partial differential equation f 1 2 . 1 2 j) for the 
effective potential we straightforwardly obtain recursive relations for even- and odd-power 
contractions W m ,k and W m +x/2,k'- 



(m—k) 

(m + h + 1) W m , k = 8 ^^;^^ + d 2 W m , k+1 + 9 2 W m ^ k 



'm-l-5/2,k-r-l 



' { 9Wl >r " 9W m ^ 2 ^ r + dW l+1/2 , r ■ dW, 

l,r 

+dW hr ■ 9W m _i_ ljit _ r+ i + dWl +1/2 ,r ■ 9F m -!-3/2,fe-r} , (3.5) 

(m—k) 

8(m + * + 2) ff„ 1Al = a '' \2k7m™-Sl^ V + * W '" +1,2 - M + # H ''"- 1/ " 

- Y { 9W l,r ■ dW m -l-3/2,k-r + dW l+1/2 ,r • dW m -l-2,k-r 
l,r 

+ dWi+i/ 2 , r ■ dW m -l-l,k-r+l + dWl, r ■ dW m -l-l/2,k-r+l] ■ (3.6) 

The diagonal contractions can easily be calculated in a closed form as in the single-particle 
one-dimensional case, yielding 

W m ,m = 1 (x-a) 2 "^, (3.7) 

(2m + 1)! 

W m+1/2jm = 0. (3.8) 



Thus, the recursion relations ( 13. 5 p and (13. 6p can be solved up to a given order p together 
with (13. 7p and (13. 8p by using a similar procedure as before. Here we give the solution up to 
order p = 4, whi ch g eneralizes the previously given solution for the simple case P = d = 1, 
obtained in Ref. [33|] . For m = we only have the naive p — 1 term, i.e. 

W ,o = V J (3.9) 

while m = 1 yields the first non-trivial even-power terms 

W ltl = \{5t-dfV, (3.10) 
o 

W lfl = ^d 2 V, (3.11) 

which are sufficient to construct p = 2 effective action. The next term we calculate is the 
odd-power contraction for m = 1, i.e. 

W 3/ 2 fl = l^-0)V, (3.12) 

which contains the explicit time derivative of V. For m = 2 we obtain the next order of 
even-power terms: 

W 2 , 2 = ^(x.fl^V, (3.13) 

W 2 ,i = T^(*" a ) 292y ' ( 3 - 14 ) 

W 2 o = — 7 + — d A V - — (dVf . (3.15) 
*° 24 240 24 v ; K J 

These terms, together with the previously calculated ones, are sufficient to construct level 
p = 3 effective action. In order to be able to complete p = 4 effective action derivation, we 
still need to calculate odd-power contractions corresponding to m = 2 

W s/2 ,i = ^{5c-d) s V, (3.16) 

1^5/2,0 = ^(x-a^y, (3.i7) 

(3.18) 



9 



as well as even-power contractions for m = 3, 



^3,3 = ^(x-9)V, (3.19) 

m ' 2 = sko (x ' 9)492 v ' (3 ' 20) 

-^(^•(x-^W, (3.21) 

^3 = JL 9 6 V + — <9 2 1> - — (didV) ■ (didV) - — (dV) ■ d 2 dV . (3.22) 
' 6720 480 360 v ; v ' 120 v ; y ' 

This concludes the calculation of level p = 4 effective action for a general many-body quan- 
tum system. As before, the obtained results automatically reduce to the already known 



effective actions for the time-independent potentials 37| if we set all time- derivatives of 



the potential to zero. Furthermore, the many-body results f l3.9D -( l3~22l) reduce for the spe 



cial case P — d — 1 to the previous time- dependent results 



331 ] . The outlined procedure 



continues in the same way for higher levels p. We have automatized this procedure and 



implemented it in our SPEEDUP code |5l| using the Mathematica software package 50| for 
symbolic calculus. 

In order to numerically verify the developed expressions for the case of mult i- component 
quantum systems, we will calculate transition amplitudes of a set of time-dependent har- 
monic oscillators, 

U(q,t) = f^(tH 2 , (3.23) 
i=i 

which represents the archetypical model for many physical phenomena. This exactly solv- 
able model allows us to compare analytical expressions obtained from recursive relations and 
to verify that using level p effective action leads to values of transition amplitudes which 
are correct up to order £ p+ 1 - pd / 2 . As we can see in Fig. [1] for a system of P = 2,4,6 
time-dependent oscillators, the respective scaling is perfect. The middle and bottom plots 
illustrate another important feature of the short-time expansion for multi-component sys- 
tems: as the number of components of the system Pd increases, the exponent p + 1 — Pd/2 
may become zero or negative for a given effective action level p. This leads to the peculiar 
behavior observed in the middle and bottom plots for small values of p, with the deviation of 
the amplitude being constant (P = 4, p = 1 in the middle plot, P = 6, p = 2 in the bottom 
plot) or even increasing (P — 6, p — 1 in the bottom plot) when e is decreased. Thus, 

10 




FIG. 1. Deviations of diagonal transition amplitudes |A^4 p (a = l,t a = 0; b = l,tj = e)| from 
the exact values as a function of propagation time e for a multi-component system (|3,23p of time- 
dependent harmonic oscillators: (top left) P = 2 oscillators, with oj\{t) = 1+ | sin 2 2t, w|(i) = 
1 + \ cos 2t; (top right) P = 4 oscillators, with w 2 (t), w 2 (t), wf(i) = 2 + cos 5i, = 4 + sin 2 At; 

(bottom) P = 6 oscillators, with w 2 (t), u^(t), w 2 (t), w|(t), w 2 (t) = 2 + sin 2 t, w 2 (t) = 4 + 2cos3t. 
Each plot gives results for transition amplitudes calculated using the effective action levels p = 
1,2,..., 16, from top to bottom. 

this prevents the calculation of the transition amplitude with high accuracy, which is, in 
principle, expected to be possible by decreasing e. As we see, to enable such high-accuracy 
calculations, one has to use an effective action with sufficiently high level p. The important 
contribution of the presented approach lies in the fact that it offers a systematic formalism 
for deriving such higher order expressions for a general quantum multi-component system. 

As in the case of single-particle one-dimensional systems, computational complexity of 
higher-order effective actions might significantly increase for higher values of p (typically 

11 



exponentially, see e.g. Fig. 5 in Ref. 33]). This increase strongly depends on the type of 
potential, but we expect to obtain significant computational benefits at least for moderate 
values of p. The optimal value of the level p to be used in practical applications can always 
be found through a small-scale numerical complexity study, by measuring relative increase 
in CPU times as a function of p. Using this data, as well as the known scaling of Monte 
Carlo errors and CPU times on the time step and number of samples, optimal value of p is 
obtained by minimizing the CPU time for a required precision of transition amplitudes. 
At the end of this section, we emphasize that the obtained discretized effective actions 



can be used for solving a p 
the exact diagonalization 



ethora of non-equilibrium many-body quantum problems within 

m 



471 ] or Path Integral Monte Carlo approach, including the 



continuous-space worm algorithm 52j . For instance, in typical experimental setups with ul- 
tracold quantum gases harmonic or anharmonic confining potentials are generically switched 
on and off, thus generating natural non-equilibrium situations. As so far mainly quenched 
potentials have been considered, it would certainly be rewarding to study in a systematic 
way how the time scale, upon which a potential is switched off, influences the observed 
time-of-flight absorption pictures. Another upcoming research field is the investigation of 
the phenomenon of parametric resonance in Bose-Einstein condensates. A first experiment, 
where the s-wave scattering length of 7 Li atoms has been modulated periodically with the 
help of a Feshbach resonance, has recently been performed [48]. In order to understand 



49] 



the observed resonance spectrum both analytical methods from nonlinear dynamics 
and numerical methods as the presented fast converging path-integral approach have to be 
combined. 



IV. VELOCITY-INDEPENDENT PART OF THE EFFECTIVE POTENTIAL 

Now we turn our attention to the special case of the velocity independent part of the 
effective potential. It determines the diagonal amplitudes, for which the discretized velocity 
x is equal to zero. The efficient and precise calculation of diagonal transition amplitudes is 
essential for many quantum statistical problems, since it provides a direct method to obtain 
partition functions, and can be used to calculate energy spectra, density profiles, and other 
relevant physical quantities. Therefore, we will derive a new set of recursive relations for 
the coefficients which determine the velocity-independent part of the effective potential. For 

12 



simplicity, we will present the derivation for one-dimensional systems, where the effective 
potential for x = can be written as 



W (x-e,r) = W(x,0;e,r) = J2 



c mfi (x,r)e m . (4.1) 



m=0 



turn out to be much simpler than the 



33|. 



Such recursive relations for coefficients c m = c m $ wil 
full set of recursions obtained in the previous paper 

In order to derive equations determining the coefficients c m , we have to perform the 
limit x — > in the partial differential equation (I2.12p for the effective potential W. This is 
nontrivial, since the equation contains derivatives with respect to x. Therefore, we have to re- 
examine both Schrodinger equations f 1 2 . X j) and (12.111) . and express them using the variables 
x, x, e, and r. After a change of variables, we get the following system of equations for the 
transition amplitude: 

1 1 

A(x,x;e,r) = 0, (4.2) 



Or 



^dd + V + - V- 



A(x,x;e } r) = 0. (4.3) 

If we take the derivative with respect to x of the second equation, use it to express the term 
dd 2 A, and insert it into the derivative of the first equation with respect to x, we obtain the 
partial differential equation 

o 4 4 4 2 2 

(4.4) 

in which it is easier to perform the required x — > limit. In the terms that do not contain 
derivatives with respect to x, we can just set x = and replace the transition amplitude A 
with A = exp (-eW )/V2^e. In the remaining terms the limit has to be performed more 
carefully. The terms containing combinations V+± VL and their derivatives are the simplest, 
and we obtain 



(V+ + V. 
d{V+ + 



m=0 

oo 

m=0 



(2m) 

V 



(2m)!2 2m ' 

e 2m+l y 

(2m + l)!2 2m+1 



B(V+ + V. 



m=0 



( 2m ), 
£ 2m y l 

(2m)!2 2m ' 



(4.5) 
(4.6) 
(4.7) 
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where the prime in the last expression denotes the derivative with respect to x. For the 
terms dA and d T dA we have to explicitly use the full double power expansion for the effective 
potential 

oo m 

W(x, x; e, t) = E r) e m ~ k x 2k + c m+1/2 , k (x, r) e m ~ k x 2k+1 ] , (4.8) 

m=0 fc=0 

perform the differentiation and take the limit afterwards. This yields the results 



dA 
d T BA 



-eA Cm+i/2 e m , (4.9) 

m=0 

oo 

-eA em (cm+i/2 - c m+ i /2 eW ) , (4.10) 



x^O 



x^O 

m=0 



where dots now represent derivatives with respect to the time argument r of the coefficient 
Cm+i/2 = c m+ i/2,o and the effective potential Wo- As we see, the odd-power coefficients c m +i/2 
cannot be eliminated altogether, although we are considering the diagonal amplitudes, for 
which we have x = 0. This is due to the derivatives with respect to x. From this we can 
deduce that we will need again two recursion relations to determine all needed coefficients, 
despite the fact that in the end we will use only the even-power ones. Therefore, we will use 
Eq. (14. 3 p to derive the second recursive relation for the coefficients. In order to do so, we 
still need to calculate the term ddA in the considered limit x — > 0: 



ddA 



A, J2 £m ( C -+V2 £ X ~ 4 + i/2 e) • (4.11) 



x^O 

m=0 



Inserting all calculated x — > terms into equations (I4.3P and (I4.4p . as well as using the 
expansion (14. ip . we finally obtain the following coupled system of recursive relations for the 
coefficients c m and c m+ i/ 2 : 

(m) i ^ ^ (20 

(2m + 1) 4 = 11(0, m) + - C-x + lfm-1/2 ~ 2 ^ ^yy^y 4-a-i 



{21+1) 

V 



I ^ >' I I 

3 . 1 . 1 

jJ2 C 'l C m-l-2 "nE C «+l/2 Cm-l-2 + J ^ cj c' r 4_j_ r _ 3 , (4.12) 



^ / j i ni—i—i 2 / j •t-v' n« »— «. • ^ 

I I Lr 



(m+1) 

2 4+1/2 = -2 n(0, m) (m+1) , 2 m+i + dm + 2 E Q +V2 4-2-1 • ( 4 -!3) 
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The above recursive relations are solved in a similar way as before. In order to obtain the 
level p diagonal effective action Wo, we need to take into account the terms in the expansion 
with m = 0, 1, . . . ,p — 1. The recursions for c' m and c- m+1 , 2 are easily solved starting from 
m = up to a desired level p—1. Although we get in this way only their first derivatives with 
respect to x, the coefficients themselves can be calculated by direct symbolic integration, and 
all solutions can be expressed in a closed form. The explicit calculation of the coefficients 
to high orders yields the same results we have obtained in the previous section. To order 
p = 1, we only have the trivial equation 

J = V, (4.14) 

that gives the well-known boundary condition cq = V. To order p = 2 we have 

c ' 1/2 = -2V + 2c = , (4.15) 
34 = \^ + \c 1/2 - 2Vc' + 2c' c = , (4.16) 

which yields c\m = and C\ = eg/ 12 = V"/12. To order p = 3 we first calculate the 
odd-power coefficient, 

c' 3/2 = 2c 1 + c 1/2 c = ^c / , (4.17) 

which leads to the result C3/2 = c /6 = V'/6. The even-power coefficient c% is obtained from 

5c' 2 = l -V' + + ^c 3/2 - 2Vc[ + ^Vc 1/2 + 4c ci + 2 C ; c - ^c 1/2 c - ^c' eg 

5V 5V /2 V , N 

— + , 4.18 

24 48 24 J ' v y 

which finally yields 

c 2 = — + - — . (4.19) 

24 240 24 

These results coincide with the results obtained in our earlier paper 



33j from the general 



recursion relations. We similarly proceed to calculate higher-order coefficients. This is 



50| , and we have 



easily automated in symbolic calculus software packages like Mathematica 
implemented the derived recursions as a part of our SPEEDUP code. 

The main advantage of this approach is that we have been able to derive recursion re- 
lations involving only the lower- level coefficients c m = c TO) o and c m+ i/ 2 = c m+ i/ 2)0 , thus not 
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FIG. 2. Deviations of diagonal transition amplitudes |A p (l,0; l,e) — .4(1,0; l,e)| as a function of 
propagation time e for: (left) time-dependent harmonic oscillator (|4.20p . calculated analytically 
for p = 1, 2, 3, 4, . . . , 20 from top to bottom; (right) forced harmonic oscillator (|4.21|) with uj = 1 
and Q = 2, calculated analytically for p = 1,2, 4, 6, 8, 10, 12, 14, 16, 18, 20 from top to bottom. The 
dashed lines on both graphs are proportional to e 15 and e 20,5 , and demonstrate the perfect scaling 
of the corresponding level p = 1 and p = 20 results. 

requiring the calculation of all even- and odd-power coefficients, which are needed for the 
general case. At the end of this section, we note that for time- independent potential V 
the recursive relation ( 14.12ft reduces to the previously known result [37j, while the second 
recursion (I4.13P yields the expected result c m+ i/ 2 = 0. Note that our recursive approach 
thus allows to calculate higher orders of the seminal Wigner expansion 531 ] . 

Fig. [2] illustrates practical advantages of using velocity-independent effective actions for 
the numerical calculation of diagonal transition amplitudes. The plot on the left gives the 
deviations of diagonal amplitudes calculated with different levels p of the effective potential 
Wq for the Grosche-rescaled 54j harmonic oscillator, 



V, 



G,HO 



[x,t) 



cu 2 x 2 



2(1 + t 2 ) 2 

while the plot on the right gives the analogous results for the forced harmonic oscillator 
1 



(4.20) 



Vfro(x, t) = -uj x —a; sin fit, 



(4.21) 



where Q denotes the frequency of the external driving field. Both models are exactly solvable, 
and the obtained e-scaling to exceedingly high orders p demonstrates the correctness of the 
analytically derived results. 
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V. REAL-TIME FORMALISM 



The presented approach has so far been developed within the imaginary-time framework, 
which is useful in many practical applications. However, in order to study the more relevant 
real-time dynamics of quantum systems, we have to switch back to the real-time formalism. 
One possibility would be to make all calculations in imaginary time, and then to try to 
perform an inverse Wick rotation, which might be difficult due to the oscillatory nature of 
the integrand in calculating the real-time propagator. Another, much more straight-forward 
possibility is to derive a new set of recursive relations within the real-time formalism. In 
this section we will briefly outline such a procedure. 

Reverting the imaginary-time formalism into a real-time one is achieved by replacing 
the variable t representing the imaginary time with it^ in all expressions, where now £r 
represents the real time. This includes also the replacement of the time-interval e with ien, 
and the time- midpoint r with its real-time counterpart ztr. Instead of (12.81) . the short-time 
transition amplitude is now expressed as 

A(a, t a , b, t b ) = - , 1 e **fcW*) (51) 

and the real-time version of the ideal effective action is defined as 

2 

S*(x, x; e R , r R ) = — x 2 - s R W(x., x; e R , r R ) , (5.2) 

which represents the counterpart of Eq. (12. 9p . Following the same procedure outlined in 
Sec. (Til we arrive at the real-time counterpart of Eq. (I2.12p for the effective potential, 

W+^-dW+ed e W- l -ed 2 W- l -ed 2 W-\e 2 {dW) 2 -\E 2 {dW) 2 = -(V++VJ), (5.3) 
8 8 8 8 2 

where the subscript R is dropped for simplicity. Further derivation of real-time recursion 
relations is a straight-forward task. For brevity, we will not give the explicit form of the 
recursion relations, but their Mathematica implementation is available from the SPEEDUP 



code web page [51 ]. 



We illustrate the applicability of this formalism for studying the real-time dynamics 
within the space-discretized approach 0, If we discretize the continuous space and 
replace it with a grid defined by a discretization step A, all quantities are only defined on 
a discrete set of coordinates q n = nA, where n G Z Pd is a vector of Pd integer numbers. 

17 



Matrix elements of the evolution operator, 

U nm (t a -)• t b ) = (q m \U(t a -)• t 6 )|q n ) , (5.4) 

represent real-time transition amplitudes A nm (t a ,tb) = A(q n ,t a ; q m ,t b ), which can be cal- 
culated using the real-time effective action approach. If the initial state of the system t a ) 
is represented by a vector if)(t a ) whose elements are if) n {t a ) — ^(Onj^o) — (QuIV'j On its 
dynamics can be calculated by a simple matrix multiplication ip(tb) = U(t a — >• £&) -tp(t a ), i.e. 

^n(*b) = ^J^nm(*a h) 1p m (t a ) ■ (5.5) 
m 

Therefore, since the matrix elements of the evolution operator can be accurately calculated 
using the effective action approach, we are able to study real-time dynamics of the system 
starting from any desired initial state. 

Note that, although we rely on the short-time expansion of transition amplitudes, we are 
not limited to study only a short-time evolution, since the above matrix multiplication can 
be repeatedly performed. For any given propagation time T, we can divide the evolution to 
iV sub-intervals of length e = T/N, which now correspond to short-time evolution matrix 
elements U nm (e). In addition to this, using the higher-order effective actions makes it 
possible to perform high-accuracy calculations of U nm , and, correspondingly, to eliminate 
the associated numerical errors for all practical purposes. 

To demonstrate this, we show in Fig. [3] the time evolution of a harmonic oscillator V(q) = 
| ui 2 q 2 calculated using the described method with effective action levels p = 1, 4 and p = 20. 
As we can see, using the propagation interval e — 0.1, the naive p — 1 action can be used 
only for short propagations times, while we are able to reproduce very accurately the long- 
time evolution of the system with higher p levels. In order to further quantitatively assess 
numerical errors of the obtained results, we use the following integral measure, 

/ poo \ 1/2 

\m)-M t )\\=[ mq,t)-4> p (q,t)\ 2 dq) , (5.6) 



where ip{q,t) represents the exact time evolution of the wave function and ip p (q,t) is the 
approximate time evolution calculated using level p effective action. The semi-log plot in 
Fig. H] gives the p-dependence of the above-defined integral measure and demonstrates that 
it obeys the expected power law, i.e. e p+1 ' 2 in this case, leading to a much smaller error 
when the propagation interval is reduced from e — 1 to £ = 0.1. This graph also shows 
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\ oj 2 q 2 calculated using 



the space-discretized method |46l. l47| and the effective action approach with: (left) p = 1, 4 and 
(right) p = 20. Both graphs display the time dependence \ip(q = 0,t)\ of the absolute value of the 
wave function at q = 0, and solid line represents the exact solution. The harmonic frequency was 
uj = 1, the time-interval for propagation was e = 0.1, and the initial state was set to the ground 
state of the harmonic oscillator with oj = 2. 

that errors due to the repeated matrix multiplication accumulate linearly with the number 
of time steps. 

The study of errors presented in Fig. [4] is very instrumental in optimizing numerical 
parameters in practical applications. If we compare errors, which correspond to the same 
total evolution time t = 1 and are calculated for a propagation in one time-step e = t and 
in N = 10 steps e = t/N, we can see that decreasing e substantially reduces errors. This 
is easily understood, since errors are proportional to £P+ l ~ Pd / 2 and, therefore, introducing 
N time steps is expected to reduce errors by a factor of N p+1 ~ Pd ^ 2 . However, the fact that 
the matrix multiplication will have to be repeated iV times introduces an additional factor 
of N, thus leading to the total reduction factor of N p ~ Pd l 2 . 

As a final example, we calculate the time evolution of the time-dependent harmonic 
oscillator with the potential 



with the frequency u(t) = 1 + ^ t for p = 1 and p = 6. Fig. |5] displays the time evolution of 
the absolute value of the wave function at q = with the propagation interval e = 0.1, and 



V(q,t) = ±uj 2 (t) q 2 



(5.7) 
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FIG. 4. Integral measure (|5.6p for numerically calculated time evolution of the harmonic oscillator 
as a function of the effective action level p for different values of the propagation time t. The 
parameters are the same as in Fig. [3l 

the initial state set to 



1 _I 2_i_Ij 



7T 



1/4 



(5.8) 



As expected, the naive p = 1 effective action can only be used for very short propagation 
times, while p = 2 action gives accurate results for longer propagation times T < 15. A 
moderate level p = 6 effective action represents an even further substantial improvement 
and can be used to accurately study much longer propagation times, as can be seen from 
the inset in Fig. [5j 

We emphasize that the presented approach might be especially of interest for Path Integral 
Monte Carlo (PIMC) calculations of the dynamics of quantum systems, in conjunction e.g. 



with the multilevel blocking method 



55] . Dynamical PIMC calculations are quite difficult 



due to the dynamical sign problem, and the availability of highly accurate propagators allows 
use of relatively small number of Trotter slices, thus substantially improving the numerical 
convergence. 



20 
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FIG. 5. Time evolution of the time-dependent harmonic oscillator (|5.7[) calculated using the space- 



discretized method 



46, 



471 ] and the effective action approach with p = 1, p = 2, and p = 6. The 



graph displays the time dependence \ip(q = 0,t)\ of the absolute value of the wave function at q = 0. 
The time-dependent harmonic frequency is given by u(i) = 1 + tq t, time-interval for propagation 
was e = 0.1, and the initial state was set according to Eq. (|5.8p . 



VI. CONCLUSIONS 



In this paper we have presented a significant extension of the approach introduced in 



the preceding paper [33|, which has established a recursive procedure for calculating the 
short-time transition amplitudes for one-dimensional quantum systems in time-dependent 
potentials. This approach is generalized here to non-relativistic many-body quantum sys- 
tems with many degrees of freedom. In parallel to the approach for time-independent po- 
tentials |37J, we have introduced an ideal effective potential for time-dependent systems, 
derived the appropriate equation using the forward and backward Schrodinger equation for 
the transition amplitude, and set up an efficient system of recursive relations, which can be 
analytically solved to high orders in the short propagation time. Furthermore, we have im- 
plemented a symbolic calculation scheme for higher-order effective actions in the SPEEDUP 
code |5l|. The analytically derived results are verified by studying several models and a list 
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of possible applications of the presented method to relevant many-body quantum systems 
has been briefly outlined. 

In addition to this, we have also studied velocity-independent part of the effective action, 
which is relevant for calculating the diagonal amplitudes and partition functions. We have 
obtained a new, simpler set of recursion relations, which determine the diagonal effective 
action, and have numerically verified that it yields the correct systematic increase in the 
convergence of diagonal amplitudes for several models. 

Finally, we have also looked at how the developed formalism can be transformed from 
its original, imaginary-time setup to the real-time one. We have derived the real-time 
counterparts of equations for the effective potential and applied the higher-order real-time 
effective actions to a numerical study of the time evolution of several models using the space- 
discretized approach 0, . We have demonstrated that the presented approach can be 
successfully used both in the real-time and in the imaginary-time formalism, and that in 
both cases we obtain analytically derived improved convergence of numerically calculated 
transition amplitudes and other quantities. We point out that the presented approach 
can contribute to improving Path Integral Monte Carlo calculations of the dynamics of 
quantum systems in conjunction e.g. with the multilevel blocking method, since higher- 
order propagators enable use of fewer number of Trotter slices. 
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